Finite Size Scaling of Mutual Information: A Scalable Simulation 
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We develop a quantum Monte Carlo procedure to compute the Renyi mutual information of an 
interacting quantum many-body system at non-zero temperature. Performing simulations on a 
spin-1/2 XXZ model, we observe that for a subregion of fixed size embedded in a system of size L, 
the mutual information converges at large L to a limiting function which displays non-monotonic 
temperature behavior corresponding to the onset of correlations. For a region of size L/2 embedded 
in a system of size L, the mutual information divided by L converges to a limiting function of 
temperature, with apparently nontrivial corrections near critical points. 



Correlation functions have been an essential tool in 
determining properties of quantum systems. However, in 
exotic phases, traditional correlation functions may prove 
insufficient to capture hidden properties. As a result, the 
Renyi entanglement entropies have attracted intense in- 
terest recently for their abilities to deliver scaling terms 
with universal numbers, e.g. at phase transitions [ij and 
in exotic topologically ordered phases [2], regardless of 
basis or choice of observable. The generalized Renyi en- 
tanglement entropies, defined by 



Sn{A) = 



1 



1 



■ln[Tr(p:^)] 



(1) 



(where pA is the reduced-density matrix of subregion A) , 
have been successfully measured recently using T — Q 
projector QMC in the valence-bond basis, by employing 
the expectation value of a "Swap" operator [3[. However, 
to study more general physics using Sn{A), one would like 
to develop a measure of entanglement entropy applicable 
for models away from SU(2) symmetry and at T > 0. 

At non-zero temperature, the mutual information (MI) 
provides the appropriate analogue of the entanglement 
entropy, measuring information between one part of the 
system and another [J|. In this paper, we develop a 
procedure for estimating the Renyi entropy in a finite- 
temperature QMC simulation, and use it to calculate a 
Renyi MI, defined in analogy to the von Neumann MI: 



In{A:B) = 5„(A) + 5„(S) - S^{AB). 



(2) 



Here, Sn{AB) is the n-th Renyi entropy for the whole sys- 
tem and Sn{A),Sn{B) are the Renyi entropies for sub- 
systems A and B respectively. At T > 0, In{A:B) is 
expected to pick up both classical and quantum correla- 
tions, while at T = 0, /„(A:B) = 25„(A) = 2Sn{B). 

Stochastic Series Expansion (SSE) quantum Monte 
Carlo ISftZl is a state-of-the-art finite-T method, able to 
efficiently simulate systems of particles with a large vari- 
ety of U(l), SU(2) or SU(N) Hamiltonians. In this paper, 
we use SSE to calculate the Renyi MI in the anisotropic 
XXZ model on a two-dimensional square lattice, with 



Hamiltonian, 
H 
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We verify our QMC procedure explicitly for small sys- 
tem sizes by showing that it reproduces l2{A:B) calcu- 
lated with exact diagonalization (ED). Then, by scaling 
the results to large system sizes with QMC, we demon- 
strate that l2{A:B) effectively captures correlations at 
both T > and T = 0. We observe non-monotonic be- 
havior in the MI, including a peak at high temperature, 
likely associated with a freezing out of classical correla- 
tions, followed by a rise at low temperature, due to the 
development of long-range quantum entanglement. Fi- 
nally, at the Ising transition for A > 1, we observe novel 
scaling behavior of the MI, where a clear crossing appears 
at Tc for different system sizes, possibly signaling a new 
type of universal scaling for this quantity. 

Finite- T QMC measurement of 5*2 . — The replica trick 
18[ for calculating the generalized Renyi entropies has 
been used extensively in field-theory, both analytically 
[l[ as well as in Monte Carlo simulations of lattice gauge 
theories [9|, llO[ . The essential formulation is to consider 
the d + 1 dimensional simulation cell with a modified 
topology - that of an n-shceted Riemann surface, with 
the region A being periodic in n/3 (in imaginary time), 
while the complement region B is periodic in /3. Restrict- 
ing to n = 2 (see Fig.[T|), the partition function Z[A, 2, T] 
of the modified system is related to the trace of the re- 
duced density matrix squared. 



Trpi 



Z[A,2,T] 
Z[TY ' 



(4) 



where Z{T) is the partition function of the "regular" un- 
modified system at temperature T. 

One immediately sees that the Renyi entropy will be 
accessible to any finite-T quantum Monte Carlo proce- 
dure that is based on a direct implementation of the mod- 
ified partition function Z[A,2,T]. At a given tempera- 
ture, 5*2 can be calculated for example through a ther- 
modynamic integration from high temperatures, analo- 
gous to the calculation of free energy in classical systems. 




FIG. 1: (color online) A schematic representation of the sim- 
ulation cell for Z[A, 2, T] in finite-T quantum Monte Carlo. 
The one-dimensional basis is represented by the horizontal 
lines, and the vertical direction is the expansion, or imaginary 
time dimension. World-lines (particle paths) are periodic in 
region A over the time length of size 2/3 (blue). In the region 
B (shaded), world lines are periodic in time /3. 



In this case, integration of two separate simulations are 
needed; one for Z[A, 2, T], and another for Z{T): 

S2{T) = -\-aZ[A,2,T] + 2\nZ{T) 



SA{fi^O)+ / {E)A,pdp 



(5) 



+ 2Sq{P = 0) - 2 / (S)o,^d/3. 

"'0 



Here, the subscript A^ /3 refers to the energy estima- 
tor calculated at a given temperature T = 1//3 in the 
Z[A,2,T] simulation, and the subscript 0,/3 refers to 
Z{T). Note that the entropy SqW = 0) = iVln(2) for an 
A^-site spin 1/2 system. For a system vifith A^^ sites in 
subregion A, Sa{I3 = 0)^[Na + 2{N- Na)] ln(2). Sim- 
ulation results used in this paper access 5'2(T) through 
Eq. ([5]); other methods such as extended ensemble tech- 
niques can also be used. 

SSE Implementation. — The above procedure can be 
implemented in any flavor of finite-T QMC based on the 
partition function, but we describe now a specific im- 
plementation in SSE QMC, referring the reader to the 
literature for the particular details and notation of the 
method Q. The SSE is based on a power-series expan- 
sion of the partition function. The essential modification 
of the usual SSE formulation is simply the direct im- 
plementation of the modified boundary configuration for 
Z[A,2,T], as described in the above illustration (Fig.[T|). 
The 5^ basis state is decomposed into spins occurring in 
region A and its complement B: \a) ~ \a'^)\a-^). The 
SSE simulation cell is composed of this basis, and a list 
of operators that propagate |a) through imaginary time. 



To simulate Z[A,2,T], two separate operator fists are 
employed, composed of ni and n2 operators, which are 
each allowed to fluctuate independently. Analogous to 
the typical SSE formulation, two separate cut-off vari- 
ables are used, Mi > ni and M2 > n2- Three crucial 
changes are thus needed in the the SSE sampling of the 
operator list and basis, as described below. 

First, note that the probability to insert or remove a 
bond operator Hi, in the diagonal update is modified: 



add 



l3Nb{ar\Hb\ar 
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(M^ - n^) 



P„ 



{M, 
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(6) 
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where b is the lattice position of the bond operator (which 
can occur on Nt, nearest- neighbor bonds in Eq. ([3])), 
and T is the "time" position in the expansion direction 
(t e [1, A/i + M2].). The index 7 has values 1 and 2 
for Z[A,2,T]. The two operators lists are sampled inde- 
pendently in the diagonal update [6J, thus their respec- 
tive size (rii or 712) is allowed to fiuctuate independently. 
Since the inverse temperature /3 occurs in each proba- 
bility above, the total number of operators ni + rii is 
on average related to the total expansion direction 2/3 
(which also suggests the straightforward extension to the 
n-sheeted Riemann surface, where 7 = 1 • • • 77,). 

The second significant modification to the SSE sam- 
pling occurs in the loop update. In fact, if one follows 
the standard two-step procedure in the loop update algo- 
rithm - first constructing a linked list of connected "ver- 
tices" and second, performing loop updates in this linked 
list y - then the modification reduces to a straight- 
forward reconnection of the topology of the linked list. 
Following intuition, the links between vertices are modi- 
fied such that, if a spin in the vertex occurs in region ^, 
connections are made as usual for the 2'(T) system with 
a periodic boundary at 2/3 (|a^) in Fig. [T]). For spins 
in region B. vertices containing operators in the list n\ 
arc connected only to themselves across the boundary 
|q;q ) at /3 (or, for vertices containing the operators in 
772, at |af ) - Fig. [1]). Remarkably, once this modified 
topology is implemented in the linked- list, the loop up- 
dates themselves are carried out as usual. In particular, 
any vertex transition weights which are used in the usual 
simulation .^(T) can be used without modification in the 
double-sheeted simulation Z\A, 2, Tl , including determin- 
istic, heat-bath, or directed-loops [6[. 

The final modification of note is the measurement of 
the expectation value of the energy in Z[A,2,T], which 
must be re-derived along the lines of Ref. [7| . Then, one 
finds that the energy is related to the total number of 
operators in both lists. 
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FIG. 3: (color online) Mutual information of the XY model 
for a square subregion A of constant size 3x3. 



FIG. 2: (color online) Measurements of Ii{A:B) (from ED) 
and hiA-.B) (from ED and SSE) for values of A = 0, 1,4 in 
Eq. Q. In the lower plot (A = 4) the exact T = values of 
Ii{A:B) (higher) and I2{A:B) (lower) are shown. The jump 
close to T = is caused by the near degeneracy of the lowest 
two energy levels of this system. The inset shows the Ii(A:B) 
curve on a linear scale, excluding the T = value, for A = 4. 



where the additional factor of 2 in the denominator can 
be modified to n for the generalized Renyi entropies. 

Simulation results for MI. — In Fig. [2l we show the Ml 
for a system of size 4-by-4 with Hamiltonian Eq. ([S]). 
where region A is a square of size 3-by-3. The precise 
correspondence between SSE and ED data for I2 clearly 
verifies the above method. Note that with ED, we are 
also able to compute /i . Two features are observed which 
deserve discussion below. 

The first feature is the low temperature behavior of the 
MI for A > 1. In the infinite system size limit, for any 
ll| . the system spontaneously breaks a discrete 



A > 1 

Ising symmetry at low temperature. As a result, for finite 
system sizes, the two lowest eigenvalues of the Hamilto- 
nian, Eo,Ei are almost degenerate {Ei — £^0 ~ 0.00009 
for this system), with a gap £'2 — Ei to the rest of the 
spectrum. For Ei — Eq <^ T <^ E2 — Ei the density 
matrix pab for the whole system is close to a mixture of 
the two lowest states V'o , V"! ■ 



PAB-{l/2){\i^0){H + \i'l){H 



(9) 



and so Si{AB) is very close to ln(2). However, since 
SiiAB) = at T = 0, the entropy Si{AB) jumps by 
« ln(2) over a very short interval of temperatures close 
to T = 0. In contrast, the entropies Si{A) and Si{B) are 
smooth at low temperatures. Hence, for A > 1 the MI 
Si{A) + Si{B) - Si{AB) jumps by approximately ln(2) 
at low temperatures due to breaking of the discrete sym- 
metry - as shown in Fig. [2] for A = 4. Thus, extracting 
the T = entanglement entropy from T > data is very 
difficult, requiring the study of very low temperatures 



for A > 1. For A < 1, the splitting between low- lying 
eigenvalues is much larger, making it easier to extract 
the T = entanglement entropy, but the presence of 
low-lying "tower of states" modes still suggests that a 
projector approach [3| should be preferred to obtain zero 
temperature data. 

A second interesting feature observed is a maximum in 
I2 at finite-T for A = 0, 4. For A = 4, this is the global 
maximum, while for A = 0, the maximum is followed 
by a further increase at low temperature. Physically, we 
interpret these maxima as follows: at very high temper- 
ature, the system has no correlations and hence no MI. 
As the temperature is initially lowered, the MI increases. 
However, eventually the fluctuations in the system start 
to freeze out, leading to a decrease in the MI. Finally, at 
low temperature, long-range quantum entanglement be- 
gins to appear, leading again to an increase in the MI. 
Note that one can also see a small maximum in the MI 
of /i for A = 4 in the inset at the bottom of Fig. [5] 

Scaling Behavior. — In Fig. [3] we plot the MI for a sub- 
region of size 3-by-3 embedded in a larger L-by-L system 
for A = 0. As L increases, the curves approach a limiting 
function of temperature - observed also for the other A 
values studied in this paper. We also consider the case 
where region A has size L/2-hy-L/2 and the system has 
size L-hy-L. In Fig. 21 we plot l2{A:B)/L as a function 
of L in the case of A = 4. The inset shows l2{A:B)/L 
over a wide temperature range, while the main figure is a 
close-up of data near T = T^ ~ 2.25 {Tc was determined 
roughly by separate measurements of the specific heat in 
Z{T)). Dividing by L shows the area law behavior, with 
the curves in the inset approaching a limiting function 
of temperature. This is an important feature of the MI: 
the difference of entropies S2{A) + S2{B) — S2{AB) can- 
cels terms proportional to volume, leaving a contribution 
proportional to the surface area of A. 

As one can see from the inset, the curves cross at two 
different temperatures. The high temperature crossing, 
at T2 ~ 4.45, can be understood simply on theoretical 




FIG. 4: (color online) Mutual information scaled by the linear 
system size L, for the XXZ model with A = 4 and a square 
subregion A of size L/2 x L/2. The inset shows two crossing 
points for I2{A : B)/L. The main figure shows detail of the 
lower-T crossing. 



grounds. For T ^ Tc, the correlations in the system 
are short-ranged, so we expect that l2{A:B) = Lf{T) + 
g{T) + C'(cxp(— L)), where the term g{T) is due to a 
correction to the MI from the four corners of region A^ 
while Lf{T) represents a contribution from the edges of 
region A (indeed, we verified that this linear fit is very 
accurate for T near T2). Then, the crossing corresponds 
to a change in sign of g{T). 

The lower temperature crossing is more remarkable, 
since it occurs at T w Tc- Our explanation of the crossing 
at high-T relies on microscopic details (the particular sign 
oi g{T)), so there does not at first seem to be any reason 
for the low-T crossing to occur near the critical point. 
One might have dismissed this as a coincidence, but we 
also observe for A = 2 the low-T crossing again occurs at 
T K, Tc- To explain this behavior, we conjecture that the 
MI as a function of temperature near Tc has the form 



h{A:B) = Lf{T) + L'^kiiT - T,y L) + g{T) 



(10) 



where < k < 1 (possibly this term depends logarithmi- 
cally on L instead). In this case, near the critical point 
the term L'^ dominates the term g{T) and the crossing 
is determined not by microscopic properties such as g{T) 
but by universal properties. 

Interestingly, the phenomenon of crossing points is also 
observed in the other models that we study. For the 
Heisenberg case (A = 1), the upper crossing occurs at 
T2 sa 0.65. There exists crossings also below the peak in 
l2{A:B)/L, (around T w 0.3 for the L studied), however 
the crossing point appears to decrease in temperature 
with increasing system size. It would be interesting to 
determine if Ti — > as A — ;■ 1, following the trend in 
Tc. For the XY model, T2 ss 0.85, while the lower tem- 
perature crossing occurs around Ti « 0.41 - a slightly 



higher value than the Kosterlitz-Thoulcss temperature 
Tkt = 0.34303(8) ^. A thorough study of the drift of 
this crossing point as a function of T, to see if it corre- 
sponds to TxT in the L ^ 00 limit would be interesting. 

Discussion. — We have presented a simulation method 
to measure mutual information at T > 0, based on the 
Renyi entropy Sn, in interacting quantum many-body 
systems. The procedure uses quantum Monte Carlo 
based on a modified partition function, with efficient lin- 
ear scaling in the spatial (A'') and temporal (n/3) simu- 
lation cell sizes. Using SSE QMC, we examine the finite 
size scaling behavior of the n — 2 mutual information 
in the spin 1/2 XXZ model with various anisotropics, 
and observe interesting non-monotonic behavior related 
to the onset of classical and quantum correlations. The 
method has allowed us to extract novel behavior at crit- 
ical points in the model, raising interesting theoretical 
questions regarding the universal scaling properties of 
the mutual information, including the possibility of a new 
critical exponent k. The generality of the idea for mea- 
suring mutual information, coupled with the versatility of 
finite-T QMC methods based on the partition function, 
should allow us to study a plethora of interesting ques- 
tions in the immediate future, including scaling behavior 
in a quantum critical fan. 
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